Slater's exchange parameters a for Analytic and Variational Xa calculations 
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Recently, we formulated a fully analytical and variational implementation of a subset of density 
functional theory using Gaussian basis sets to express orbital and the one-body effective poten- 
tial. The implementation, called the Slater-Roothaan (SR) method, is an extension of Slater's Xa 
method, which allows arbitrary scaling of the exchange potential around each type of atom in a 
heteroatomic system. The scaling parameter is Slater's exchange parameter, a, which can be de- 
termined for each type of atom by choosing various criteria depending on the nature of problem 
undertaken. Here, we determine these scaling parameters for atoms H through CI by constraining 
some physical quantity obtained from the self-consistent solution of the SR method to be equal to 
its exact value. Thus, the sets of a values that reproduce the exact atomic energies have been de- 
termined for four different combinations of basis sets. A similar set of a values that is independent 
of basis set is obtained from numerical calculation. These sets of a parameters are subsequently 
used in the SR method to compute atomization energies of the G2 set of molecules. The mean 
absolute error in atomization energies is about 17 kcal/mol, and is smaller than that of Hartree- 
Fock theory (78 kcal/mol) and the local density approximation (40 kcal/mol) but, larger than that 
of a typical generalized gradient approximation (~ 8 kcal/mol). A second set of a values is de- 
termined by matching the highest occupied eigenvalue of the SR method to the negative of the 
first ionization potential. Finally, the possibility of obtaining a values from the exact atomization 
energy of homonuclear diatomic molecules is explored. We find that the molecular a values show 
much larger deviation than what is observed for the atomic a values. The a values obtained for 
atoms in combination with analytic SR method allow elemental properties to be extrapolated to 
heterogeneous molecules. In general, the sets of different a values might be useful for calculations 
of different properties using the analytic and variational SR method. 

PACS numbers: 

Keywords: density functional theory, exchange potential, Slater's Xct 



The Hohenberg-Kohn-Sham (HKS) formulation of 
density functional theory is by far the choice for to- 
day's electronic structure calculationiSiiii Prior to the 
HKS formulation of the density functional theory (DFT), 
Slater formulated its basis in an attempt to reduce 
the computational complexity of the Hartree-Fock (HF) 
method. 3 Unlike the HF method, this Xa method has an 
exchange potential v x that is local, and proportional to 
the one-third power of the electron density p. It is given 
by 



Vx [p] 



3/3\i/3 1/3 



where a is called Slater's exchange parameter and was 
originally equal to one. Gaspar— and Kohn-ShamS ob- 
tained similar expression on other rigorous grounds. 
They variationally minimized the total energy functional 
and determined a to be 2/3, in contrast to Slater's value. 
The difference in the two values is rooted in the aver- 
aging process employed in the simplification of the HF 
exchange potential. Slater, taking a cue from Dirac's 
earlier work^ averaged the exchange potential over the 
entire Fermi sphere of radius kf = (3n 2 p(r)) 1 ' 3 while in 
the Gaspar's and Kohn-Sham procedure the value of ex- 
change potential at k = kf (r) is used. It was suggested 
that the a parameter in Slater's method could be treated 



as an adjustable parameter^ The set of a parameters for 
atoms that make self-consistent Xa energy match the 
HF energy was determined by Schwarz. 8 The a values 
he obtained typically range from 0.77 for light atoms to 
0.69 for heavy atoms. He also noted that different atomic 
configurations lead to only slight changes in the a values, 
suggesting that the Xa method can also be applied to the 
molecules or solids. Several other ways to determine the 
a parameters have been put f rthi24i2iiiii2ii&±2ii£ii& 

Early implementations of the Xa method for molecules 
with atom-dependent a parameters employed the muf- 
fin tin (MT) approximationpiiiS*i2iSliSi In this scheme, 
atoms or ions are enclosed by the atomic spheres in which 
the potential is approximated to be spherically symmet- 
ric while, in the interstitial region, the potential is con- 
stant. The effective one body potential in the MT imple- 
mentation is discontinuous at the surface of MT sphere 
and therefore the energy was mathematically undefined. 
This method, however, had an advantage over all the 
other quantum-chemical methods which was that the 
molecules dissociated correctly. At infinite inter-atomic 
distances a sum of atomic energies could be reproduced. 
When the MT model was discarded due to difficulties in 
geometry optimization, a single value of a ( = 0.7) had 
to be chosen for all types of atoms S There are, however, 
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some attempts to obtain good thermo-chemistry from the 
Xa calculations by correcting or improving the uniform 
a calculations in a secondary calculation. 23 ' 24 25 ' 26 

Electronic structure calculations using traditional 
quantum chemical methods such as the HF theory, and 
beyond, generally employ basis sets to expand molec- 
ular orbitals. In these methods, calculations of ma- 
trix elements and other quantities of interests are an- 
alytic. On the other hand, almost all implementa- 
tions of density functional models, until recently, re- 
quired use of numerical grids to compute contribu- 
tions to matrix elements from the exchange-correlation 
terms. Cook and coworkersSSiS&S 9 * 3 ^ 3 ^ successfully 
demonstrated fully analytic implementation of density 
functional model. His variational implementation used 
Gaussian basis sets to express molecular orbitals, the 
one-body effective potential^ and the uniform a Xa ex- 
change potentialiSSiS The advantages of analytic calcu- 
lations are obvious. The calculations are fast and accu- 
rate to machine precision (within the limitations of the 
model and basis set). Round-off error, which grows as 
the square root of number of points and other numeri- 
cal problems are eliminated^ Consequently, smooth po- 
tential energy surfaces are obtained. 28,29 Unfortunately, 
modern sophisticated exchange-correlation functionals 
are too complex to allow fully analytic solution at this 
time. Analytic DFT, at this stage, is restricted to Slater's 
exchange-type functionals. We have recently proposed an 
algorithm that permits fully analytic solutions and also 
allows for atom dependent values of Slater's exchange 
parameter a for heteroatomic systems^ This method, 
called the Slater-Roothaan (SR) method, is based on 
robust and variational fitting and requires four sets of 
Gaussian bases. It can have the advantage of MT Xa 
that molecules dissociate correctly as interatomic dis- 
tances tend to infinity. 3 * Intuitively, it is natural to 
expect that when molecules dissociate the constituent 
atoms would have exact atomic energies at infinite sepa- 
ration. This important physical requirement is satisfied 
within the SR method if a parameters that reproduce 
exact atomic energies are used in molecular calculations, 
although one might want to use other sets of a's for other 
molecular properties. 

In this work, we report different sets of a parameters 
that could be used in the SR method. The first such 
a set is determined by requiring that the self-consistent 
atomic SR atomic energy be equal to the exact atomic 
energies. The a values in this set will be hereafter re- 
ferred to as atomic alpha values. We have recently used 
some of the a values from this set to calculate the total 
energies for the G2 sei£^ of molecules and found that 
these a values give remarkably good total energies 
Here, we use them to compute the atomization ener- 
gies for the G2 set of molecules and compare them with 
other models. Further, we also obtain two additional 
sets of a based on other criteria for possible use in the 
SR method. The second set of atomic a values is de- 
termined to provide negative of eigenvalue of the high- 



est occupied molecular orbital (HOMO), chomo, that 
matches with the exact value of the first ionization po- 
tential are determined. 37 i 38 i 39 i 40 i 41 i 42 The third set of a is 
determined to provide the exact atomization energy for 
the selected homonuclear diatomic molecules of the first 
and second row of the periodic table. All calculations 
are repeated for four different combinations of analytic 
Gaussian basis sets. 

It is apparent that the analytic Slater-Roothaan 
method is empirical in nature as it uses adjustable 
a values. It, however, differs from the semi- 
empirical^ 3 ^^ methods which use a minimal basis set, 
or empiricai 4 ^^^ 8 ^^ models which do not compute elec- 
tronic structure, and the minimal-basis tight binding 
methods^ reported in the literature. It is a variant of 
density functional models that is also computationally 
very efficient because it is fully analytic and requires 
no numerical integration. We have successfully used it 
for studying the heterogeneous systems such as boron 
and aluminum nitride nanotubes containing up to 200 
atomsi£i*££ Furthermore, it also provides the flexibility of 
being tuned, through the a parameters, according to the 
need of problem. The sets of a parameters reported in 
this work will be useful starting points in this regard. In 
the following section we outline the analytic SR method 
used in this work. In section II, we describe computa- 
tional details which are followed by results and discussion 
in section III. 



I. THEORETICAL METHOD 

The total electronic energy in the DFT for an TV- 
electron system is a functional of electronic density p and 
is given by 

JV 

eKS \p\ = E < > + E *AphPi\ (!) 

i 

where, the first term contains the kinetic energy operator 
and the nuclear attractive potential due to the M nuclei, 

The total electron density is expressed in terms of the 
Kohn-Sham orbitals </>j i<T as 

p(r)=p- t + p l , (3) 

with 

i 

where rij )CT is the occupation number for the <pi i(T or- 
bital. The second term in Eq. Q represents the classical 
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Coulombic interaction energy of electrons 

p(r)p(P) 



E„ 



(p\\p) 



d 3 rd 3 r'. 



(5) 



This energy is approximated by expressing the charge 
density as a fit to a set of Gaussian functions, 



(6) 



where, ~p(f) is the fitted density, di is the expansion coef- 
ficient of the charge density Gaussian basis function Gi. 
The elimination of the first order error in the total energy 
due to the fit leads toSl 



E ee = 2(p\\p)-(p\\p). 



(7) 



The expansion coefficients {di} are determined by varia- 
tionally minimizing this energy with respect to {di}. The 
last term F, xc in Eq. is the exchange energy, 



Exc[p-[,Pi\ 



6\V3 



d 6 r 



pf(r)+pf(r) . (8) 



The form of above functional allows analytic calculations 
with the Gaussian basis to be performed. For this pur- 
pose the one-third and two-third powers of the electron 
density are also expanded in Gaussian basis sets: 



p*(r) fa pa = ^e t E t {f) 



(9) 



pt(0«p* =J2m(f). (10) 



Here, {Ei} and {Fi} are independent Gaussian basis 
functions, while and /j are expansion coefficients. The 
exchange energy is then given bjiSiSLifflii^iiS 3 . 



Exc — Ca 



^(pph - ^(p* p* ph + \{p' 3 ph 

1/3 



(ii) 



where C' a = — 9a(^j Thus, using the four LCGO 

basis sets (one for orbital expansion and three fitting ba- 
sis sets) the total energy is calculated analytically. The 
LCAO orbital coefficients and the vectors d, e, and f are 
found by (constrained if desired) variation. 

A. Slater-Roothaan method 

The expression for the total electronic energy in the 
Slater-Roothaan method has the following form: 33 

E SR = Y / <<t>i\h\<l>i>+2(p\\p)-(p\\p) 

i 

[J(flB3?> - 



<r=1,l 



+ 3(^6 S 3 fi) 



(12) 



TABLE I: The optimal a values for the different basis sets 
that yield the exact atomic energies in the highest symme- 
try for which the solutions have integral occupation numbers. 
The numerical values are for the spherically symmetric atoms 
and are obtained by the numerical atomic structure code. The 
exact atomic energies given in the last column are from RefloTl 
The basis sets are I: 6-311G**/RI-J, II: 6-311G**/A2, III: 
DGDZVP2/RI-J, and IV: DGDZVP2/A2. 

Basis I Basis II Basis III Basis IV Numerical E (a.u.) 



H 

Li 

Be 

B 

C 

N 

O 

F 

Na 



Si 
P 

s 

CI 



0.77739 
0.79169 
0.79574 
0.78675 
0.77677 
0.76747 
0.76500 
0.76066 
0.75204 
Mg 0.74994 
Al 0.74822 
0.74539 
0.74324 
0.74262 
0.74197 



0.77739 
0.79169 
0.79574 
0.78668 
0.77672 
0.76747 
0.76495 
0.76067 
0.75204 
0.74994 
0.74819 
0.74540 
0.74324 
0.74260 
0.74196 



0.78124 
0.79211 
0.79614 
0.78684 
0.77670 
0.76726 
0.76454 
0.76002 
0.75287 
0.75120 
0.74872 
0.74600 
0.74397 
0.74352 
0.74273 



0.78124 
0.79211 
0.79614 
0.78677 
0.77665 
0.76726 
0.76448 
0.76001 
0.75287 
0.75120 
0.74869 
0.74602 
0.74397 
0.74350 
0.74272 



0.77679 
0.79118 
0.79526 
0.78744 
0.77657 
0.76654 
0.76454 
0.75954 
0.75110 
0.74942 
0.74797 
0.74521 
0.74309 
0.74270 
0.74183 



-0.500 
-7.478 
-14.667 
-24.654 
-37.845 
-54.590 
-75.067 
-99.731 
-162.260 
-200.060 
-242.370 
-289.370 
-341.270 
-398.140 
-460.200 



Here, the partitioned 3/4 power of the exchange energy 
density, 



,9B 



(f) = C*5>(*)a(j)£> 



(13) 



where C x — C a /a and -D?-(r) is the diagonal part of the 
spin density matrix, and the function, 



a(i) = a 



3/8 



(14) 



contains the c^, the a in the Xa , for the atom on which 
the atomic orbital i is the centered. 



II. COMPUTATIONAL DETAILS 

The a values are calculated for the atoms that oc- 
cur in standard set of (Becke's)^ 56 molecules. These 
atoms are H, Be through F, and Na through CI. Mag- 
nesium and aluminum atoms do not occur in the Becke's 
set but have been included in the present work. The 
set of 56 molecules actually consist of 54 molecules with 
two molecules in two different electronic states. The SR 
method requires four basis sets, of which one is the orbital 
basis. The remaining three basis sets are required for fit- 
ting of the the Kohn-Sham potential (See Ref. |3l| for 
more details). For the orbital basis sets, our choices are 
the triple-C (TZ) 6-311G** basisS^ and the DGaus^l 
valence double-C (DZ) basis se& called DZVP2. The 
most reliable and thus best s-type fitting bases are those 
that are scaled from the s- part of the orbital basisiSl 

o 1 

The scaling factors are two for the density, # for pa 
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and | for ~pi . These scaled bases are used for all s-type 
fitting bases. A complete package of basis sets has been 
optimized 58 for use with DGauss. 57 In addition to the va- 
lence double-^ orbital basis, called DGDZVP2 herein, we 
use the pd part of the (4,3;4,3) for Be-F (A2) charge den- 
sity fitting basis. Ahlrichs' group has generated a RI-J 
basis for fitting the charge density of a valence triple-^ or- 
bital basis set used in the Turbomole program. 59 These 
fitting bases are used in combination with the 6-3 11G** 
and DGDZVP2 basis sets. The a are obtained for the 
different combination of above basis sets. The molecule 
geometries were optimized using the Broyden-Fletcher- 
Goldfarb-Shanno (BFGS) algorithmi 60 i 61 i 62 i 63 i 64 Forces 
on ions are rapidly computed non-recursively using the 
4-j generalized Gaunt coefficientsiS 5 ^& The calculations 
are spin-polarized for open-shell systems. 



III. RESULTS AND DISCUSSION 

The optimal a values that give the exact atomic ener- 
gies can be obtained by the Newton-Raphson procedure 
that finds zeros of the function f(a) — E SR {a)—E exact — 
0, where E and E exact are the self-consistent Xa (SR) 
and the exact total energies, respectively. Finding zero of 
this function requires frequent calculation of this function 
and its derivative. Determination of derivative, /'(a), is 
straight forward. Only the exchange term in the energy 
functional E SR (a) depends explicitly on the a parame- 
ter. The derivative obtained as the exchange-energy di- 
vided by the a parameter provides sufficiently accurate 
approximation to the actual derivative. The LCAO fits 
depend weakly on a. The Newton-Raphson procedure 
was implemented using PERL scripts and the energy and 
its derivative with respect to a were obtained from the 
FORTRAN90 SR code. The sets of a values obtained 
using this procedure are given in Table in|U It is evident 
from the table that the choice of fitting basis does not 
significantly affect the a values. The a values do show 
some dependence on the orbital basis set. However, the 
changes in alpha values are small with the a values for 
the 6-311G** basis set being consistently smaller than 
those for the DGDZVP2 basis set for the same fitting 
basis set. For the same choice of fitting basis sets, the 
6-311G** orbital basis set will give lower energy than 
the DGDZVP2 basis. Therefore, the a values for the 
DGDZVP2 basis should be larger than the 6-311G** ba- 
sis to provide the more negative (or binding) potential 
required to give the energy equal to the exact energy. 
Except for H, which has no correlation energy, these a 
values are larger than the reported HF values^ which is 
expected as these values are obtained for the exact en- 
ergy (bounded from above by the HF energy). The HF 
a values show systematic monotonic decrease with the 
atomic number^ The current a values show overall de- 
crease with the atomic number except for peaks at Li and 
Be. The peak at Li is caused by the fact that H has no 
correlation energy while peak for Be is caused by electron 



TABLE II: The optimal a values for the different basis sets 
for which the eigenvalue of the highest occupied molecular 
orbital is exactly equal to the exact ionization potential (IP) . 
The experimental values of the first ionization potentials (in 
eV) are also included in the last column. The basis sets are I: 
6-311G**/RI-J, II: 6-311G**/A2, III: DGDZVP2/RI-J, and 
IV: DGDZVP2/A2. 

Basis I Basis II Basis HI Basis IV IP^ 

H 1.1901 TT901" 1.1877 1.1877 13W 
Li 1.1246 1.1246 1.1152 1.1152 5.39 
Be 1.2749 1.2749 1.2581 1.2581 9.32 
B 1.1227 1.1220 1.0993 1.0989 8.30 
C 1.1006 1.1000 1.0775 1.0770 11.26 
N 1.0832 1.0832 1.0629 1.0629 14.52 
O 1.2360 1.2345 1.1958 1.1946 13.61 
F 1.1696 1.1694 1.1501 1.1498 17.42 
Na 1.1181 1.1181 1.1221 1.1221 5.14 

Mg 1.2256 1.2256 1.2311 1.2311 7.64 
Al 1.1113 1.1112 1.1074 1.1067 5.98 
Si 1.0864 1.0836 1.0875 1.0844 8.15 
P 1.0501 1.0501 1.0662 1.0662 10.49 
S 1.1543 1.1562 1.1759 1.1770 10.35 
CI 1.1211 1.1218 1.1344 1.1347 12.96 



pairing and correlation in the outer orbital. 

The atomic a values can also be determined by differ- 
ent criteria. In the KS DFT it has been argued that the 
negative of the eigenvalue of the highest occupied orbital 
—tmax equals the first ionization potentiali2LiSiiSi2Mii22i 
Although such an interpretation of the — e m ax is not yet 
settled, it gives us another way to determine set of a pa- 
rameters. Such a such set of a's might be useful in future 
for calculations of polarizability or optical spectra. These 
a values are given in Tablc[H]for the four combinations of 
basis set. All the values presented are larger than those 
obtained from the total energy matching criteria. This 
is not surprising as the — e max for the local density type 
(a = 2/3) functionals typically underestimate the ioniza- 
tion potential by as much as 30-50%. It is known that 
this occurs due to the incorrect asymptotic behavior of 
the effective potential. The asymptotic behavior of the 
effective potential is governed by the exchange potential. 
In the present Xa method, the exchange potential decays 
exponentially in the asymptotic region, due to which the 
valence electrons experience shallower long-range poten- 
tial than they otherwise should. In the present work, we, 
however, treat a as purely a fitting parameter to get the 
right e max . These a values do not decrease left to right 
across the periodic table and are rather close to Slater's 
value of a. The removal of self-interaction of electron also 
leads to better asymptotic description^ We think that 
the self-interaction correction can be also be implemented 
in the present analytic SR method. Orbital densities are 
non-negative. 

Finally, we explore the possibility of obtaining exact 
atomization energy (AE) within the present SR method, 
for selected diatomic molecules. For this purpose we use 
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TABLE III: The optimal a values for the selected homonu- 
clear dimers that yield the exact atomization energies for the 
6-311G**/RI-J basis set. The exact formation energies and 
the bond lengths (in A) are also included in the last column. 





Ot 


Do (kcal/mol) 


Re 
present 


Re 

Expt. 7 


Ho 


1.39172 


-103.5 


0.59 


0.74 


Li2 


1.45747 


-24.0 


2.12 


2.67 


Be 2 


0.56596 
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2.46 


c 2 
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N 2 
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Si 2 
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-100.7 
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Newton- Raphson procedure to obtain the zero of the fol- 
lowing function: f(a) = E? n R ol {a) - 2E s at R om {a) -AE = 0, 
Here, E^ t is the self-consistent total energy of the op- 
timized molecules for a given a, and Ef R m is the self- 
consistent total atomic energy for a given alpha. The 
bond length of the molecule is optimized during each 
Newton- Raphson step. These calculations were per- 
formed using PERL scripts and the calculated set of a 
values is presented in the Table ITTT1 These values show 
large variation as function of atomic number. The gen- 
eral trend is that the a values for the dimers of atom on 
the left of periodic table are larger than 1.0 and smaller 
than 1.0 for the atoms on the right. For the exact val- 
ues of the atomization energies the bond lengths of these 
dimers show the general trend that the atoms on the left 
side of the periodic table are larger and those on the right 
are smaller. The exception are the Be2 and N2 dimers, 
perhaps due to their weak and very strong bonding, re- 
spectively. Another noteworthy observation is that these 
dimers with exception of Be2 are still bound for vanish- 
ingly small value of a . This is interesting as no molecules 
are bound in the Thomas-Fermi models It seems the 
bonding therefore occurs due to the exact treatment of 
kinetic energy in the present Xa method. Here, it should 
be noted that the 6-311G** basis sets used in this work 
may not provide satisfactory description for very small a 
values as the electrons then will experience a very shallow 
potential and will be rather delocalized leading thereby 
to the artifacts like positive eigenvalues for the outer- 
most electrons. For this reason the a for the F2 dimer 
should be used cautiously. Further, these a values when 
viewed relative to the Kohn-Sham's value of 2/3 provide 
some hints on how the corrections to the Kohn-Sham 
exchange should energetically contribute for accurate de- 
termination of atomization energies. The molecules with 
molecular a values larger than the Kohn-Sham value of 
2/3 will be underbound in the exchange only KS scheme 
and those with a lower will be overbound. Any univer- 



TABLE IV: The atomization energies Do (kcal/mol) for the 
56 set of molecules calculated using the two sets of a values 
given in Table I for the two basis sets that reproduce the exact 
atomic energies. The two basis sets are I: 6311G**/RIJ, II: 
DZVP2/A2. The last column contains the exact values. 
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TABLE V: The mean and mean absolute error (kcal/mol) in the calculated atomization energy of 56 molecules (Cf. Table 
llVfl relative to their experimental values for different basis sets. The SR (EA) are present calculations while SR (HF) are SR 
calculation with olhf values. The SR (HF) results are from Ref. l33t The mean absolute error for the standard Hartree-Fock 
theory and the local density approximation is 78 and 40 kcal/molf— *— respectively. OB: Orbital basis, FB: Fitting basis 
(pdfg-type). 



OB 


FB 


SR{EA) 


SR(EA) 


SR(HF) 


SR(HF) 






mean 


absolute 


mean 


absolute 


6-311G** 


RI-J 


-1.9 


17.3 


-5.1 


16.4 


DGDZVP2 


A2 


-1.6 


16.2 


-4.8 


16.1 



sal correction to the KS exchange functional should be 
such that it simultaneously reduces overbinding in some 
molecules and underbinding in others. 

The atomization energies for the set of 54 molecules 
in 56 the electronic states are presented in the Table llVl 
These are calculated by the SR (EA) method using the 
set of a values given in Table I. The calculations are per- 
formed for the two sets of basis sets. These molecules 
in the dissociation limit will give the exact atomized en- 
ergies. It is evident from the Table ITvl that the Xcv un- 
derbinds some molecules while it overbinds others. This 
trend is in contrast with the local density approxima- 
tion and the HF theory. The former uniformly overbinds 
while the latter uniformly undcrbinds. A similar trend 
has been found when the HF values of a are usedi^ In 
general, the error in atomization energy is smaller when 
the molecules consist of atoms on the opposite sides of 
the periodic table, for example, HC1, H2S, NaCl, HO etc. 
The best agreement is observed for the NH3 with the 
DGDZVP2 basis and OH with all basis. The mean and 
mean absolute errors in the atomization energies of these 
molecules calculated with respect to their experimental 
values is given in Table [V] These errors are somewhat 
larger than those obtained when the HF a values are 
used. We note here that the HF value of a, cchf, are 
those that when used in the Xa model reproduce the HF 
exchange energies for atoms. 

The a values obtained from the e max matching crite- 
ria can be used to obtain the first ionization energy of 
a molecule from its highest occupied eigenvalue. Here, 
we demonstrate their use for calculations of the first ion- 
ization potential of N2 and CO molecules. For the N2 
molecule, using the a values from Tabic [H]for the basis-I 
through basis-IV, the ionization potential is 15.38,15.34, 
15.40, and 15.20 eV, respectively. These values are in 
good agreement with experimental value of 15.58 eVii 
For the CO molecule, the ionization potential is 14.83, 
14.83, 14.53, and 14.53 eV for the basis sets I through IV. 
These values are also in very good agreement with CO's 
experimental ionization energy (14.01 eV)H We expect 
that use of these a values to give better eigen-spectra, 
which could be used in the the polarizability calculations 
by sum over states method^ The atomic or Hartree- 
Fock a values could be used to optimize the geometry 
of molecule and then to calculate the matrix elements 
required in the sum over states expression. The eigen- 



values required in the sum over states can be obtained 
by performing one more self-consistent calculations us- 
ing the a values from Table ^ The use of improved 
eigen-spectrum in the sum over states method has been 
found to give good estimates of polarizability J2> Calcu- 
lations for the set of 54 molecules using molecular a val- 
ues (Cf. Table IIII|) are not performed because the SR 
method as implemented now^ can not handle wide vari- 
ation in the molecular a values. Moreover for the rea- 
sons mentioned earlier these a values (particularly for 
the molecules containing F) necessitates bigger orbital 
basis sets. The trend in the bond lengths in Table IIIII 
suggests that for these dimers better values of a may 
be obtained by minimizing the deviation of both the dis- 
sociation energy and the bond length, from their exper- 
imental values. There are several other possibilities for 
determination of a which may be useful in improving the 
performance of the SR method. As mentioned earlier 
the differences in the Slater's and Gaspar-Kohn-Sham 
value of a are the result of different averaging process 
employed in the derivation of exchange potential. This 
suggests the possibility that different a values for the de- 
scription of the core (say a c ) and valence (a v ) electrons 
may provide better description of the exchange correla- 
tions. The additional parameter in this case will provide 
greater flexibility in the fitting procedure. Alternatively, 
the Xa exchange functional can be augmented with suit- 
able functional forms that allows the analytic solution 
of the problem. The local functional form by Liu and 
coworkers that consists of the sum of 0, 1, 1/3, and 2/3 
powers of electron density and originates from the adia- 
batic connection formulation appears to be particularly 
suitable for this purpose^*^ We are currently exploring 
these possibilities. 

To summarize three sets of a values are obtained 
on the different criteria for the use in the analytic SR 
method. The first set of a values is determined by equat- 
ing the self-consistent total energy of atoms to the "ex- 
act" atomic energy. Two other sets are determined us- 
ing different constraints such as equating the negative 
of eigenvalue of the highest occupied orbital to the first 
ionization potential and reproduction of exact atomiza- 
tion energies for the diatomic molecules. The exam- 
ination of the performance of the SR method for the 
atomization energies of 56 molecules with the first set 
of a values gives the mean absolute error to be about 
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TABLE VI: Mean absolute error (MAE) in calculated to- 
tal energies of G2 set of 56 molecules for different models. 
The errors are in kcal/mol and are at optimized geome- 
tries in respective model. Ml: SR(EA)/6-311G**/RI-J, M2: 
SR(EA)/DGDZVP2/A2, LDA: local density approximation 
(See Ref. |3^ for more details) 



Model 


MAE 


Ml 


17.3 


M2 


16.2 


LDA 


532 


PBE-GGA 


101 


B3LYP-hybrid GGA 


15 



17 kcal/mol. This value is far larger than the gen- 
erally accepted chemical accuracy of 2 kcal/mol. The 
tabulation of the MAE for the G2 set and extended 
G2 set for more sophisticated functionals is given in 
RefFl The MAE is within 4-9 eV for the parameteriza- 
tions at the level of the generalized-gradient approxima- 
tion (GGA), 78 i 79 i 80 i 81 The meta-GGAS or hybrid GGA 
functionals 54-83 perform even better with MAE of about 3 
eV. The MAE is 37-40 kcal/mol for the local spin density 
approximation^ The SR method's performance is inter- 
mediate between the local spin density approximation 
and the GGA. We note that we have not optimized the 
method to give accurate estimate of any particular prop- 
erty. The a parameters could be chosen to minimize the 
MAE in atomization energies like the other density func- 
tional models, including the GGA, do. Such an optimiza- 
tion process would necessarily lower the error in Table 
IVl Our goal in this work is to examine simple schemes 
for extrapolating elemental properties to heterogeneous 
molecules, without any additional massaging of the re- 
sults. For example, the set of a values from Table [I] can 
be used to extrapolate atomic energies to molecular ener- 
gies by SR method. Indeed, such an extrapolation gives 
remarkably accurate total energies (see, Table E3, that 



are comparable to or better than those obtained by some 
popular, sophisticated pure and hybrid density functional 
models*^ By construction, using the a values of this Ta- 
ble, the MAE in total energies is same as that in atomiza- 
tion energy. This is in contrast to many popular density 
functional models which seem to give better atomization 
energies due to cancellations of errors in total energies of 
atoms and molecules. Also, the SR method is unique in 
that molecules dissociate correctly in the separated atom 
limit. It is therefore not unreasonable to use a values 
adjusted for atoms. The present work shows that fitting 
perhaps any molecular property can be done quantum 
mechanically through density functional theory. If the 
calculations are analytic and variational, then small ba- 
sis sets can be used to generate unique, stable and reliable 
energies using minimal computer time. Thus one can en- 
vision embedding quantum-mechanical calculations, with 
full geometry optimization where appropriate, to opti- 
mize very sophisticated quantum-mechanical calculations 
of molecular properties over the G2 and larger sets of 
molecules using PERL scripts to control the optimiza- 
tion on a single processor or to farm out independent 
sub-optimizations on multiple processors. Our toolbox 
of functionals that can be treated analytically contains 
more 85 than the cube-root functional used in this work. 
Thus one can expect that as more experience is gained 
in analytic DFT we can better approximate the best ex- 
change and correlation functionals that currently require 
limited-precision numerical integration. The process of 
driving MAEs lower and lower through better and bet- 
ter analytic functionals need never end, short of perfect 
agreement. 
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